Emergence of a novel sublineage, MYMBD21 under SA-2018 lineage of Foot-and-Mouth Disease Virus serotype O in Bangladesh

Foot-and-Mouth Disease (FMD) hinders the growth of the livestock industry in endemic countries like Bangladesh. The management and prevention of FMD are severely impacted by the high mutation rate and subsequent frequent generation of newer genotypes of the causative agent, Foot-and-Mouth Disease Virus (FMDV). The current study was conducted in nine districts of Bangladesh during 2019–21 to characterize the circulating FMDV strains based on the VP1 sequence analysis, the major antigenic recognition site providing serotype specificity and high variability of FMDV. This study detected the first emergence of the SA-2018 lineage in Bangladesh along with the predominance of Ind-2001e (or Ind-2001BD1) sublineage of ME-SA topotype under serotype O during 2019–21. The mutational spectrum, evolutionary divergence analysis and multidimensional plotting confirmed the isolates collected from Mymensingh districts, designated as MYMBD21 as a novel sublineage under the SA-2018 lineage. Analysis of the amino acid sequence revealed several changes in the G-H loop, B-C loop and C-terminal region of VP1, revealing a 12–13% divergence from the existing vaccine strains and a 95% VP1 protein homology, with most of the mutations potentially considerable as vaccine escape mutations, evidenced by three-dimensional structural analysis. This is the first report on the emergence of the SA-2018 lineage of ME-SA topotype of FMDV serotype O in Bangladesh, as well as a possible mutational trend towards the emergence of a distinct sublineage under SA-2018 lineage, which calls for in-depth genome-wide analysis and monitoring of the FMD situation in the country to implement a strategic vaccination and effective FMD control program.


Scientific Reports
| (2023) 13:9817 | https://doi.org/10.1038/s41598-023-36830-w www.nature.com/scientificreports/ were endemic and were not detected outside the respective countries 16,17 . Ind-2001 lineage is diversified into five sublineages: Ind-2001a, b, c, d, and e 3,6,18 . Recently, a distinct novel lineage, SA-2018 under ME-SA topotype was reported from India in 2018 19 . The high mutation rate of the viral genome, lack of proofreading activity of the RNA polymerase and the quasi-species nature of the virus are responsible for the emergence of genetically and antigenically diverse strains of FMDV 19,20 .
In Bangladesh, Ind-2001d, Ind-2001BD1 and Ind-2001BD2 under the Ind-2001 lineage of ME-SA topotype were found to circulate 3 . Later, Ind-2001BD1 was designated as Ind2001e by the World Reference Laboratory for Foot-and-Mouth Disease (WRLFMD) 6,19 . The Ind-2001e or Ind-2001BD1 sublineage was predominant over the other sublineages 3 . Ind-2001BD2 sublineage was confirmed as a novel sublineage in a previous study and circulation of this particular sublineage was not found after 2013. Few uncategorized isolates (BD_BAU_ML1_2013; BD_BAU_ML2_2013; BD_SI_5_2013; O/BAN/BLRI/450.2/2018) were reported during 2013-18 21 . BD_BAU_ ML1_2013 and BD_BAU_ML2_2013 isolates were reported as PanAsia-2 lineage 22 as the isolates shared 98-99% identity with some Indian isolates reported in 2011 but these sequences were not genetically characterized 23 . Therefore, no clear information was provided about the lineage of these two isolates.
The existence of immunogenically diverse strains of FMDV without cross-protective immunity contributes to the devastating outbreaks each year in Bangladesh. FMD outbreaks appear as a threat to the economy of Bangladesh. To prevent the disease realistically, a well-planned control strategy must be designed that requires constant surveillance of the circulating FMDV strains and measuring the effectiveness of the existing vaccines against the circulating local strains. This study focused on molecular typing of the FMDV serotype O isolates collected during 2019-21 and reported the infiltration of a new lineage from the neighboring country as well as a mutational trend towards the emergence of a novel sublineage under the recently introduced lineage in Bangladesh. Antigenic heterogeneity of the newly emerged strain against existing vaccine strains was also documented in this study. A Table S1, Table S2). The sequences from the Mymensingh district, named MYMBD21 in this study showed an evolutionary relationship to SA-2018 isolates, reported from India in 2018 (Fig. 2). MYMBD21 formed a completely distinct clade from the Indian isolates of the SA-2018 lineage with 5-6% divergence ( Supplementary Fig. S1, Fig. S2, Fig. S3). MYMBD21 also showed a relationship with previous Bangladeshi strains which were reported during 2013 (named as PanAsia-2 lineage) and 2018 (uncategorized) 21,22 but these isolates were 8% distant from both the MYMBD21 and Indian isolates of SA-2018 ( Supplementary  Fig. S4). The genetic identity of MYMBD21 was closer to SA-2018 than previous unclassified isolates of Bangladesh (

Conformational change between MYMBD21 with SA-2018 Indian isolate and uncharacterized
Bangladeshi isolate. The superimposition of VP1 structures of MYMBD21 against SA-2018 isolate and uncharacterized Bangladeshi isolate showed some minor changes in the structure due to substitutions in crucial sites but no significant conformational changes in the G-H loop of VP1 occurred. The mutated sites were labelled on the structure. The superimposed structures are given in the Supplementary information file (Fig. S20 a, b).  In the next two years, the changes in sequence identity from 99% to 97% between the circulating strains and Indian isolates was observed which was possibly driven by a selective pressure from both the geographical and demographic factors of this country ( Fig. 2; Supplementary Table S1).

Discussion
In the phylogenetic analysis, the sequences from the Mymensingh district, designated as MYMBD21 in this investigation, revealed an evolutionary link with the 2018 India-reported SA-2018 isolates and prior uncategorized isolates from Bangladesh (BD BAU ML1 2013, BD BAU ML2 2013, BD SI 5 2013, and O/BAN/ BLRI/450.2/2018) under the same ME-SA topotype (Fig. 2). However, MYMBD21 developed a completely separate clade from the SA-2018 cluster, differing by 5-6% ( Fig. 2; Table 1; Supplementary Fig. S1, Fig. S3) and also from prior uncategorized isolates from Bangladesh with 8% divergence (Fig. S4). The genetic distance of MYMBD21 with SA-2018 (0.062) was lower than that of the Bangladeshi isolates (0.089) and any of the established lineage (Table 1; Supplementary Table S3). The level of identity and genetic distance analyses both suggested that MYMBD21 belongs to the SA-2018 lineage. The genetic distance between previously circulating uncharacterized Bangladeshi isolates and the SA-2018 lineage was 0.101, which was analogous to distances among other established lineages, indicating that these uncategorized Bangladeshi isolates was not under the SA-2018 lineage (Table 1; Supplementary Table S3). This demonstrates that SA-2018 lineage has not previously been discovered in Bangladesh. According to OIE annual country report-2020 (October-December), 13 VP1 sequences of collected samples from India showed the formation of a distinct clade, called 2018 lineage which was later named SA-2018 lineage 19 . In Bangladesh, no study reported the presence of this new lineage in circulation before 2021 until confirmed by this study.
The MYMBD21 isolates, once confirmed to be under the new lineage with the formation of a unique clade, were further analysed for its credibility to be called a new sublineage. BLAST search results and genetic distance analysis showed a considerable divergence (5-6%) from Indian isolates that suggests the possible generation of a new sublineage (Supplementary Fig. S1; Table 1). In multidimensional scaling, MYMBD21 isolates formed a separate cluster from the clusters of other lineages or sublineages, yet were found to be the closest one to SA-2018 cluster (Fig. 3). Formation of a separate evolutionary clade in phylogenetic tree, genetic distance sufficient    (D85N, L126M, R140H). Two mutations happened to be prior to the G-H loop at 85 th position where negatively charged aspartate (D) was mutated to polar uncharged asparagine (N) and at 126 th position where Leucine (L) was converted to Methionine (M) both having hydrophobic side chains. One mutation occurred at the antigenically important G-H loop (R140H) where arginine (R) was transformed to histidine (H) (Fig. 4).
VP1 amino acid sequence of MYMBD21 differed by 5% from both vaccination strains, with 11 alterations (Supplementary Fig. S8, Supplementary Fig. S9; Figs. 5 (Fig. 7a, b). Changes in the polarity and charge in the amino acids in the crucial antigenic sites might be attributed to the major conformational changes in G-H loop.
The study shows the first report on the emergence of the SA-2018 lineage in 2021 accompanied by the prevalence of Ind-2001e (or BD1) sublineage during 2019-2021 in Bangladesh. The co-circulation of multiple FMDV serotype O lineages and sublineages, however, can cause comparable cross-immunity between dominant and emergent lineages, which may be crucial for lineage turnover as in endemic areas, susceptible populations develop immunity to the dominant lineage of FMDV rather than the emerging lineage 25 . The extent of a novel lineage's dominance in the future is often difficult to predict, but given the FMDV serotype O lineage's evolutionary history, the SA-2018 lineage can be potential outbreak-initiating lineage in Bangladesh. Additionally, the mutational spectrum of the lineage is evident to reflect a strong indication for the generation of the sublineage O/ME-SA/ SA-2018/MYMBD21 in Bangladesh. The impact of the SA-2018 lineage is going to last for a longer time, as the susceptible population will encounter the rapid transmission efficacy of FMDV. It is possible that significant FMD outbreaks might be caused by the novel SA-2018 lineage around the year 2022 given the temporal trend of serotype O lineages emerging and then dominating after a time period of 3-4 years observed in the field. As a result, appropriate measures, such as constant surveillance of FMD outbreaks, appropriate control program and effective vaccination program must be put in place to limit the spread of the disease.

Conclusions
Intrusion of a new lineage SA-2018 in Bangladesh followed by its emergence as a novel sublineage along with structural heterogeneity with vaccine strains complicates the current FMD situation in Bangladesh. Given the events from the recent past, this sublineage may have the potential to become dominant in near future which requires further genome wide analyses, vaccine candidate selection and stringent control measures to be adopted as soon as possible.  Fig. S21). Sample collection protocol was in accordance with guidelines and regulations approved by the AEERC and all authors complied with the relevant guidelines. All methods are reported in accordance with ARRIVE guidelines (https:// arriv eguid elines. org). All the steps of sample handling and processing were performed in a biosafety level 2 laboratory facilities.

Sample collection.
A total of 156 tongue or foot epithelium tissue samples were collected from FMD- www.nature.com/scientificreports/ Polymerase chain reaction-based amplification of VP1 and VP1 sequencing. VP1-based PCR diagnostic assay was employed for the detection of FMDV positive tissue samples. VP1 region of cDNA was amplified using two sets of primers (VP1UF/NK61, 16F/NK61). The PCR reaction was performed using GoTaq 2 × Hot Start Colorless Master Mix (Promega, USA) with either forward primer VP1UF (5′-GTA CTA CRCSCAG TAC -3′) 21,22 or 16F (5′-GAG AAC TAC GGW GGW GAG AC-3′) 12 and the reverse primer NK61 (5′-GAC ATG TCC TCC TGC ATC TG-3′) 9 . The PCR products were run on 1.0% agarose gel with a 1 kb-DNA ladder (Promega, USA) for the visualization and detection of FMD positive amplicons. Following detection, the FMDV positive amplified PCR products were purified using the Wizard SV Gel and PCR Clean-Up System (Promega, USA). Then, PCR products were sequenced using BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, USA). The VP1 coding sequence quality was analysed in ABI Genetic Analyser (Applied Biosystems, USA). Both forward and reverse sequences were assembled into a single contig using SeqMan version 7.0 (DNASTAR, Inc., Madison, WI, USA). The assembled sequences were compared with other sequences from GenBank 26 29 with the related gene sequences from GenBank. Sequences were accessed from the database on August, 2022. Based on the alignment data, best substitution model Tamura-Nei with discrete Gamma distribution (TN93 + G) was selected based on the lowest BIC (Bayesian Information Criterion) score. Then, phylogenetic Maximum Likelihood tree based on the Tamura-Nei model 30 were constructed (bootstrap replicates 1000). A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+ G, parameter = 0.6693)). Alignment gaps, missing data, and ambiguous bases of about fewer than 5% were allowed at any position as all the sequences were not completely aligned on the full range. The percentage of nucleotide identity between sequences was retrieved from the BLAST 27 search and was also calculated using a global alignment tool based on the Needleman-Wunsch algorithm 31 in NCBI.
Genetic divergence analyses of proposed novel sublineage from established lineages or sublineages and previously reported Bangladeshi FMDV isolates were calculated based on the genetic variations to confirm the emergence of novel sublineage or lineage in Bangladesh using the Kimura 2-parameter model 32 in MEGA11 28 . The rate variation among sites was modelled with a gamma distribution (shape parameter = 1).

Multidimensional scaling (MDS).
Multidimensional scaling (MDS) was used to visualize the divergence or closeness between the proposed sublineage and other established lineages or sublineages of FMDV. For this analysis, selected reference VP1 sequences from each established lineage or sublineage of serotype O were aligned along with suspected new sublineage and pairwise genetic distance was calculated using the Kimura-2 parameter model 32 in MEGA11 28 . Using the distance matrix, MDS plotting was performed using IBM SPSS Statistics (Version 26). Each point in the diagram represented one object from the distance matrix of sequences.
Analysis of VP1 amino acid variations. The VP1 amino acid sequences of samples under study, previously reported reference sequences and vaccine strains were used to analyse mutational spectrum among them. The amino acid was translated based on the standard genetic code after codon-based alignment with MEGA11 28 and BioEdit 33 was used for the visualization of variable amino acids and the positions where the substitutions occurred.
Prediction of the 3-D structure of VP1 Region. Homology modelling of VP1 of a representative FMDV sample, two reference sequences, and two vaccine strains, was performed using the SWISS-MODEL server 34,35 . STML ID 5nem.1.A 36 was used as a template for modelling the FMDV sample, and STML ID 5ner.1.A 36 was chosen to model the rest of the structures. The quality of three-dimensional structures was validated using the Ramachandran plot 37 analysis in Molprobityv4.4 38 . The quality and the energy criteria of structures were assessed by the ProSA-web 39 using the PDB files of each model. Using PyMOL 40 , the models were visualized and superimposed for necessary structural analysis.

Data availability
Supplementary materials supporting the results of the study are available in this article as Supplementary datasheet and Supplementary information file. VP1 sequences of FMDV isolates included in this study have been submitted to the GenBank databases under accession number OP320415-OP320458; OP271696-OP271697.